Gut microbiota differences in stunted and normal-lenght children aged 36–45 months in East Nusa Tenggara, Indonesia

The role of the gut microbiota in energy metabolism of the host has been established, both in overweight/obesity, as well as in undernutrition/stunting. Dysbiosis of the gut microbiota may predispose to stunting. The aim of this study was to compare the gut microbiota composition of stunted Indonesian children and non-stunted children between 36 and 45 months from two sites on the East Nusa Tenggara (ENT) islands. Fecal samples were collected from 100 stunted children and 100 non-stunted children in Kupang and North Kodi. The gut microbiota composition was determined by sequencing amplicons of the V3-V4 region of the 16S rRNA gene. Moreover, fecal SCFA concentrations were analyzed. The microbiota composition was correlated to anthropometric parameters and fecal metabolites. The phyla Bacteroidetes (Bacteroidota; q = 0.014) and Cyanobacteria (q = 0.049) were significantly higher in stunted children. Three taxa at genus levels were consistently significantly higher in stunted children at both sampling sites, namely Lachnoclostridium, Faecalibacterium and Veillonella (q < 7 * 10−4). These and 9 other taxa positively correlated to the z-score length-for-age (zlen), while 11 taxa negatively correlated with zlen. Several taxa also correlated with sanitary parameters, some of which were also significantly different between the two groups. All three fecal SCFA concentrations (acetate, propionate and butyrate) and their total were lower in stunted children compared to non-stunted children, although not significant for butyrate, indicating lower energy-extraction by the gut microbiota. Also, since SCFA have been shown to be involved in gut barrier function, barrier integrity may be affected in the stunted children. It remains to be seen if the three taxa are involved in stunting, or are changed due to e.g. differences in diet, hygiene status, or other factors. The observed differences in this study do not agree with our previous observations in children on Java, Indonesia. There are differences in infrastructure facilities such as clean water and sanitation on ENT and Java, which may contribute to the differences observed. The role of the gut microbiota in stunting therefore requires more in depth studies. Trial registration: the trial was registered at ClinicalTrials.gov with identifier number NCT05119218.


Introduction
Stunting in children is defined as those children being too short for their age, defined as a stature with a height-for-age (or length-for-age as used by the World Health Organization [WHO]) score more than 2 standard deviations below the reference standard [1].This is mostly due to poor nutritional intake [2], but also hygiene practices, including access to clean water, together with socioeconomic influences, all contribute to the health status of these children and influences growth.Stunting has become a major health problem in the world.Because a mother, who experienced stunting in childhood, will likely bear a stunted child, this is a cyclic process [3].The linear growth retardation in these children already begins in utero and continues into infancy and early childhood [4].The 2020 data from UNICEF/WHO/ World Bank Group [5] reveals that more than one-fifth of children under five-years-old, a staggering 144 million children worldwide, were stunted due to chronic malnutrition.Of these 54.3% were found in Asia, and 39.9% in Africa.This stunting was largely irreversible after the child's second birthday.
Stunting therefore has been identified by UNICEF as a major global health priority [6].In line with that, the WHO has set a target to reduce stunting by 40% between 2010 and 2025.Moreover, improving the identification, measurement and understanding of stunting is part of recommended actions by WHO in scaling up the prevention [7].Although frequently stunting is equated with malnutrition and chronic infection, the latter mostly due to lack of hygiene, there is accumulating evidence that the link with malnutrition and infection alone may not be entirely valid, as evidenced in a population of Indonesian children [8].This study ". . .seriously questions the concept of stunting as prima facie evidence of malnutrition and chronic infection" [8].
Indonesia is considered to have a high prevalence of stunting overall.According to the global nutrition report [9], Indonesia experienced a slight decrease in stunting from 42.4% in 2001 to 36.4% in 2013.Based on current data, the cases of stunting have been decreasing further to 30.8% in 2018 [10] and to 27.7% in 2019 [9], but still amounts to 4.6 million children in 2019.Although there is a trend for a gradual decrease, the burden of stunting in Indonesia is still above the rate of stunted children in the Southeast Asia region (24.7%)[9], and depends heavily on the district within Indonesia in which these children live [11].Around 53% of the population now live in urban areas, whereas the remainder live in rural areas spread over the many islands, with large differences in hygiene practices and facilities [5].
The gut microbiota has close links to food digestion, food absorption and intestinal function.Persistent undernutrition in childhood will change the (development of the) healthy composition of the gut microbiota (referred to as dysbiosis) [12].Vice versa, gut microbiota dysbiosis has been associated with malnutrition and reduced plasma essential amino acid levels in children from two low-income countries: Malawi and India [13].There is accumulating evidence that the gut microbiota influences body weight regulation, both in overweight and obesity [14] and well as in anorexia nervosa [15].The microbial activity, in terms of metabolites produced by the gut microbiota, has been shown to play a role in weight regulation, with a particular role for the short-chain fatty acids (SCFA; primarily acetate, propionate and butyrate) [16,17].Also in stunting, an altered gut microbiota is linked to the pathophysiology of stunting, where such alteration have even been detected prior to actually observing delayed growth in children between 6 and 23 months [18].In the context of stunting, is has been observed that Enterobacteriaceae (of the Proteobacteria phylum [or Pseudomonadota phylum as it has been renamed to [19]], which is often associated with human pathogenicity) are increased in concordance with impaired digestion/absorption and localized gut inflammation [20,21].Gut microbiota differences are attributable to our own genetics, diet and age, amongst many other variables.The composition is individual and region specific, with different compositions in Asia, Europe, US and Africa [22][23][24][25], partly correlated to differences in diet for Asian children [22], although that data is obtained in children at two location in 5 countries (Bali and Yogyakarta for Indonesia) and may not be representative for the whole country, or rural areas.In the Asian Microbiome Project, particular profiles in school-aged children could be linked to intake of resistant starch (e.g. from certain rice varieties) and chicken [22].
The profile of fecal microbiota of apparently healthy and those of stunted Indonesian children aged 3-5 years needs to be explored to find out how gut microbial community structure changes with nutritional status.In our earlier study in 2019 on microbiota in stunted children on Java, Indonesia [26,27], we found that the composition was correlated to nutritional status.In stunted children, the relative abundance of Bacteroidetes was significantly lower than in normal children, while Firmicutes was significantly higher.At the genus level, 14 taxa were significantly different between groups.Overall "Prevotella 9" was the most abundant genus (average of 27%), and it was significantly lower in stunted children than in normal children (23.5% vs. 30.5%,respectively).The two districts on Java that were chosen then (Pandeglang, Banten province and Sumedang, West Java province), were amongst the regions in Indonesia with had the highest prevalence of stunting at 22% and 37.8%, respectively [8], despite the fact that both districts are close to capital cities in Indonesia.Based on the Study on the Nutritional Status of Indonesia (2021) [11] East Nusa Tenggara has 15 districts with high stunting prevalence; among others, Kupang city and Southwest Sumba (26.1% and 44%, respectively).
The aim of the study was to explore the interrelationships between the gut microbiota profile and the nutritional status of children in East Nusa Tenggara (Kupang-urban versus North Kodi-rural), and to identify key microbial taxa associated with stunting.This information can subsequently be used to correct the observed dysbiosis in these stunted children by optimal interventions to manage malnutrition, through proper dietary interventions, and access to proper clean water, personal hygiene and sanitation practices, which in turn modulate the gut microbiota.

Study design
A cross-sectional study was conducted in children aged 3-5 years old with stunting (n = 100) and those without stunting (n = 100), at two locations: Kupang and North Kodi, in the East Nusa Tenggara (ENT) Province, Indonesia from November 06, 2021 to March 08, 2022.
This was a follow-up of a previous study on Java, which also include 100 stunted and 100 non-stunted children [27].The same conditions (such as sample storage, primers, PCR cycles, sequencing protocol, bioinformatic tools, etc.) were used in this study as the one performed by the same team on Java earlier.

Collection of feces
Fecal samples were collected on site from stunted children and non-stunted children and immediately placed in a cooler with ice-packs, and shipped the same day to the lab in Jakarta on dry ice.In the lab, 0.5 g of the feces was mixed with 4.5 ml of Zymo buffer (Baseclear, Leiden, the Netherlands) and kept at room-temperature prior to extraction of DNA.Another 0.5 gram was weighed in a tube for extraction of SCFAs.

DNA isolation and sequencing of the V3-V4 region of the 16S rRNA gene
DNA isolation and sequencing of barcoded amplicons of the V3-V4 region of the 16S rRNA gene were essentially performed as described before [28] according to established protocols provided by Illumina (Illumina, Eindhoven, the Netherlands).The sequencing was carried out using the Illumina MiSeq system (San Diego, CA, USA).Sequences were converted into FASTQ files using BCL2FASTQ pipeline, primer sequences were removed and quality trimming was applied based on the Phred quality score.QIIME 2 software was used for microbial analyses [29].The sequences were classified using SILVA (version 132) as a reference 16S rRNA gene database.Alpha-and β-diversity metrics were obtained using QIIME2 as well.The Amplified Sequence Variants (ASVs) table was filtered for those taxa that occurred in more than 20% of the total number of fecal samples.Blanks, positive controls (a fecal sample that is tested with each sequence run) and a mock sample (also tested with each sequence run) were taken along and behaved as expected (i.e.: no reads, the expected composition at 10.000 reads, and the expected relative abundance at 10.000 reads, respectively).

Short chain fatty acid analyses
SCFAs were analyzed by GC-MS, essentially as described before [28].Briefly, fecal samples were mixed 1:1 (w:w) with PBS, vortexed for 2 minutes and centrifuged at 14,000×g for 10 min.One hundred and fifty μl of the supernatant were mixed with 550 μl internal standard solution, containing methanol, internal standard (2 mg/ml 2-ethyl butyric acid), and formic acid (20%).The analysis was carried out on a GC-MS (8890 GC System; Agilent Technolgies, Amstelveen, the Netherlands) equipped with a PAL3 RSI 85 autosampler (Agilent) by injecting 1 μl sample on a DB-FATWAX Ultra Inert column (30 m, 0.25 mm, 0.25 μm, Agilent).The temperature settings of the injector port, oven, flame-ionization detector and mass spectrometer detector were 250, 200, 275 and 225˚C, respectively.The flow rate over the column was 1.2 ml|/min.The MS Quantitative Analysis (Quant-My-Way) software from Agilent was used to determine concentrations of SCFA.LODs were: 0.95 mmol/L for acetate and 0.91 mmol/L for both propionate and butyrate.The lowest concentrations found in the samples were at least 2-fold higher than the LOD.

Statistical analyses
Correlations between ASVs and different categorical variables, such as stunted/non-stunted, gender, sampling site, or sanitary variables were investigated using the non-parametric Kruskal Wallis test corrected for multiple comparisons by Benjamin-Hochberg false discovery rate (FDR), by using the software package R (3.5.3) (R Core Team, http://www.R-project.org/) in RStudio.Non-parametric Spearman's rank-order correlations were obtained between ASVs and continuous variables, such as age, height and weight, and the concentrations of SCFAs present in feces.Q-values (adjusted p-values after FDR) were considered significantly different at a strict cut-off of q < 0.05.The map of Indonesia with the taxa plotted in pie-charts was made with the R-package 'maps' (https://cran.r-project.org/web/packages/maps/index.html).Beta-diversity was visualized using principal coordinate analysis (PCoA) plots in QIIME2.Permutational multivariate analysis of variance (PERMANOVA; [30]) was performed to test the significance of beta-diversity (weighted and unweighted UniFrac, Bray-Curtis dissimilarity and Jaccard similarity) between normal and stunted children in QIIME2.

Ethics approval and consent to participate
The study protocol was approved by the Ethics Committee of Research Institute of YARSI University and registered at ClinicalTrials.gov with identifier number NCT05119218.Parents or children's caregivers received oral and written information and signed a letter of consent before children included in the study.

Characteristics of the children 36-45 months of age
This was a cross-sectional study, without intervention.Fig 1 shows the CONSORT flow chart of the study.Table 1 shows the anthropometric parameters of the two groups of children segregated based on their zlen.Stunted children were on average 8 cm (9%) shorter compared to the children with normal zlen (Table 1, p = 1.8 * 10 −10 ).They were also significantly lighter in bodyweight (p = 1.5 * 10 −6 ), as well as different in BMI (p = 2.6 * 10 −6 ).The age of the two groups did not significantly differ (p = 0.62).
Fig 2 shows the distribution of the zlen of non-stunted children (green) and stunted children (red), plotted against the WHO reference dataset (grey).Stunting is defined as Z-score < -2, and severe stunting as Z-score < -3.Some of the children had a Z-score lower than -4 (n = 8), with the lowest score being -5.01.From the distribution of the children considered non-stunted, it is clear that most of these (n = 86 of 100) are to the left of the zero-line in the WHO reference data-set, and therefore on average are still shorter than the average length-forage in the WHO-dataset.This was similar to our previous study on stunted children in comparison to non-stunted children on Java, Indonesia [27].

Differences in sanitary parameters between the groups
Several sanitary parameters were significantly different between stunted and non-stunted children (Table 2).The remainder of the parameters did not reach significance.
On all hygienic parameters in Table 2, except kitchen ventilation, the stunted group scored lower, meaning worse, than the normal group.This has been observed before in Indonesia [31,32], and elsewhere [33] before, reflecting not only nutritional intake, but also hygiene factors that are responsible for stunting.

Gut microbiota differences between stunted and non-stunted children
The gut microbiota of the stunted children (S) from Kupang and North Kodi were compared to those of the non-stunted children (N) from the same region.Kupang is the capital city and the administrative center of East Nusa Tenggara.It is the largest city in East Nusa Tenggara  and an urban area with the population of approximately 443,000 people.It is categorized as a semi-arid area due to relatively low rainfall, tropical, and dry climate which also tends to be influenced by the wind.North Kodi on the other hand is a rural area consisting of mountains and hills with a long dry season.The North Kodi community lives relatively far from urban areas.
With respect to α-diversity, the Shannon-and Faith's phylogenetic diversity (PD-)indices were different between stunted children and non-stunted children (q = 0.084 and 0.0072, respectively).Evenness (q = 0.245) and observed_features (q = 0.056; which is called observe-d_OTUs in QIIME2, but measures ASVs) were not significantly different (Table 4), although the latter showed a trend to be higher in stunted children.The α-diversity was also compared between the urban and rural sampling sites, and evenness was significantly different between sampling sites (Table 4).This led us to split up the samples and compare stunted versus normal per sampling site.Like for the complete dataset, for samples collected in Kupang, the Shannonand Faith's PD-indices were different between stunted children and non-stunted children (Table 4).
However, none of the α-diversity-indices were significantly different for samples collected in North Kodi.Moreover, there were significant differences in Shannon-, Faith's PD-and observed_features-indices between samples collected from stunted children in Kupang compared to those of children in North Kodi.Observed_features was also different between samples from non-stunted children from Kupang and North Kodi (Table 4).With respect to β-diversity, except for the unweighted UniFrac between samples from stunted and non-stunted children in North Kodi, all other metrics were significant, indicating that there were not only differences between children being stunted or not, but also between collection sites and the interaction between these two (Table 4).Given the differences observed between the two sampling sites (urban and rural), we wondered whether the Kupang children might cluster closer to the samples we collected earlier in two cities on the western part of Java (Pandeglang and Sumedang, in Banten and West Java provinces, respectively) [27].The β-diversity comparisons show significant differences between these two studies (Table 4), and the Bray-Cutis dissimilarity PCoA plot shows a very clear separation between the samples from Java (Fig 3B ; open circles) and those from East Nusa Tengarra (closed circles).Similar separation between the two datasets are also observed for the other β-diversity metrics (S2 Fig) .Also, some differences in α-diversity between the two dataset are eminent (Table 4).The total number of taxa at the genus-level in the full dataset was 412.After filtering for those that were present in at least 20% of the samples, 125 genus-level taxa remained, indicating that the majority of genera was present in only a few children.Of the 125 taxa remaining, 52 (41.6%) were significantly (q < 0.05) different between sampling sites using the nonparametric Kruskal-Wallis (KW) analysis, possibly influenced by previous observed differences between rural and urban microbiota composition [35,36], observed already as early as 2001 [37].We did not investigate this further, as we were interested in the differences correlated with stunting.Taking the full dataset (all 200 children), KW analysis indicated 30 taxa to be significantly different between stunted children and non-stunted children (q < 0.05).However, after segregating the data into the rural and urban sampling sites, only 3 taxa were consistently different between stunted children and non-stunted children at both sampling sites (Fig 4 ): Lachnoclostridium (q = 3.1 * 10 −6 ), Faecalibacterium (q = 3.1 * 10 −4 ), and Veillonella (q = 6.7 * 10 −4 ).
In our previous study, Prevotella 9 was the most prevalent taxon, with on average 26.3% relative abundance (RA) [27], and it was lower in stunted children on Java (23.5% for stunted children vs. 30.5% for non-stunted children; q = 0.058).In the current study however, the RA of Prevotella 9 was much lower (on average 9%), and opposed to our previous results was higher in the stunted children than in the children with normal nutritional status (9.95% for stunted children vs. 8.04% for non-stunted children; q = 0.36).Prevotella has previously been linked with non-Western subjects [38].However, it has also been linked to ingestion of (plant) carbohydrates [39].How Prevotella relates to stunting exactly remains to be seen given the disparity in results from our two studies.
In the current study Faecalibacterium was the taxon with average highest RA (11.1%; 13.3% for stunted children vs. 8.8% for non-stunted children; q = 3.1 * 10 −4 ; Fig 4B1).Faecalibacterium was slightly higher in the children from North Kodi compared to those from Kupang (Fig 4B2 ).Interestingly, Faecalibacterium is considered to be a biomarker for health, with associations with numerous diseases and disorders, such as IBD (particularly Crohn's disease) [40] and obesity [41].It is therefore counterintuitive to find higher RA of Faecalibacterium in stunted children.
Lachnoclostridium has previously been associated with autism spectrum disorder (ASD) in children [42], and is used as a novel marker for colorectal cancer in adults [43].It also has been shown to be involved in production of secondary bile acids in critically ill children [44].Veillonella are opportunistic pathogenic bacteria.The overgrowth of Veillonella is associated with many diseases, including oral diseases, endocarditis, and ulcerative colitis, and an overabundance of Veillonella promotes intestinal inflammation [45].Veillonella has also been shown to be modulated in ASD, but not consistently, with some studiesshowing increased Veillonella and others decreased Veillonella in ASD [42].
The second highest taxon in the current study was Escherichia/Shigella.It was highly prevalent, 9.7% in stunted children and 9.3% in non-stunted children, but not significantly different between the two groups (q = 0.71).Strikingly, in children in East Nusa Tenggara (ENT) (9.5%) this taxon was approximately 7-fold higher than in children on Java (1.4%; the green color can hardly be observed in Fig 5 for the Java samples).This may be related also to observed differences in diarrheal prevalence in these provinces: diarrhea prevalence in ENT, Banten and West Java is 9.1%, 8.4%, 5.4%, respectively [11].The same data from the Indonesian ministry of health also revealed that access to proper drinking water varies between these provinces: in ENT, West Java, and Banten it is 46.6%, 67.6%, and 72%, respectively.Moreover, proper sanitation in ENT is also lower than on Java: with 62.9% in ENT versus 82.1% and 75.2% in Banten and West Java, respectively [11].Differences in culture, education, infrastructure facilities are also quite diverse between Java and Nusa Tenggara islands, and these may contribute to the gut microbiota profile of the young children.
The three taxa that were significantly different between the two groups of children in East Nusa Tengarra were not significantly different between the two groups of children on Java, and vice versa.Strikingly, Veillonella was only present in 1 child in the Java-samples.Although there was a difference in average age of the two studies (average 46.8 months within Javastudy, average 39.7 months within the ENT-study), and although age does influence gut microbiota composition as reviewed by Selma-Royo et al. [46], none of the three taxa different between stunted children and non-stunted children were correlated with age, not in the current study, nor if both datasets were combined.
Other studies have looked at the microbiota differences between stunted and normal children.In another region in Indonesia, on the northern tip of Sumatra, the gut microbiome profle of stunted children (n = 21) showed enriched genera such as Blautia, Dorea, Collinsella, Streptococcus, Clostridium sensu stricto 13, Asteroleplasma and Anaerostipes.Meanwhile, compared to the normal children (n = 21), there was a depletion of Prevotella, Lactococcus, Butyrivibrio, Muribaculaceae, Alloprevotella, Akkermansia, Enterococcus, Terrisporobacter and Turicibacter [47].In a small longitudinal analysis of the gut microbiota in persistently stunted young children in South India (n = 10 children in each group), the RA of the Bacteroidetes phylum was higher in stunted children compared to that of control children at 12 months of age but otherwise did not find significant differences in diversity indices between cases and controls.The microbiota of control children was enriched in Bifidobacterium and Lactobacillus species, whereas that of stunted children was enriched in inflammogenic taxa including those in the Desulfovibrio genus and Campylobacterales order [48].In another study on South-Indian children [49], stunted children (n = 16) showed enrichment for bacterial genera including Prevotella 7, Prevotella 9 and Sutterella.They also showed depletion of Clostridiaceae 1 family, Bifidobacterium, Intestinibacter and Fusicatenibacter genera compared to non-stunted children (n = 25).In 75 infants aged 5-12 months living in rural areas of Peru that were followed for 6 months, the RA of several bacterial taxa were found to increase over time in children that became stunted (n = 16) but not in non-stunted children (n = 49).These included Ruminococcus 1, Ruminococcus 2, Ruminococcaceae UCG 014, Lachnospiraceae, Clostridiaceae 1, and Collinsella.Providencia decreased over time in children that became stunted [50].
In the only other large study, in 335 children from rural Zimbabwe from 1-18 months of age, the authors found, using whole metagenomic sequencing, that taxonomic microbiome features are poorly predictive of child growth, however functional metagenomic features, particularly B-vitamin and nucleotide biosynthesis pathways, moderately predict both growth and growth velocity.The prevalence of stunting varied from 18-34% across study time-points [51].
A recent review describes decreased microbial diversity and changes in specific bacterial taxa in malnourished children.These alterations have been linked with compromised gut barrier function, immune dysregulation, and systemic spread of the gut microbiota, but changes in composition between different studies were not consistent [52].

Correlations between gut microbiota and sanitary parameters
Several taxa correlated with sanitary parameters (Fig 7A).Type of toilet was correlated to several taxa (Fig 7A ), but not significantly different between the two groups (Fig 7B).However, whether it was public or private, and located inside the accommodation or elsewhere were important discriminators between the two groups (Table 2 Interestingly, the three taxa that were significantly different between the two groups, Lachnoclostridium, Faecalibacterium and Veillonella, were not among the taxa that correlated with the hygienic parameter 'toilet type' ( Fig 7).Instead, Akkermansia, a taxon generally considered to be correlated with health, did show a correlation with toilet type (q = 0.0087).How this relates to stunting remains to be seen, as it was not correlated to being stunted or not (q = 0.26).Rinanda et al. [47] noticed that water contaminants such as Aeromonas, Stappiaceae, and Synechococcus were higher in stunted children compared to normal children.In our study, an overall increase in the Escherichia/Shigella taxon, usually correlated with hygiene, was observed (not specifically in stunted children), as well as a specific increase in stunted children in the phylum Cyanobacteria, of which the natural habitat is considered to be water environments.

Correlations between gut microbiota and microbial metabolites
Fecal SCFA concentrations were measured as well.Although we acknowledge that these concentrations cannot be interpreted into SCFA production by the microbiota, it is interesting to observe the correlations between the taxa and the presence of these microbial metabolites in feces.It is difficult to convert SCFA concentrations to SCFA production, because first of all, if concentrations between the groups is the same, but fecal output (which was not recorded) differs, then the absolute amount of excreted SCFA differs too (despite the same concentrations).Moreover, an increased fecal excretion can be caused by a lower intestinal absorption, and thus does not necessarily mean increased production.Despite this, all three SCFA concentrations were lower in stunted children, compared to non-stunted children, although only significant for acetate and propionate (Table 5).As a consequence, also total SCFA was also lower in stunted children than in children with normal zlen (p = 0.011).If this could be translated to lower SCFA production, there is less energy-extraction from the diet by the gut microbiota in the stunted children.Also, since SCFA have been shown to be involved in gut barrier function [14,17], this barrier may be negatively affected in the stunted children.Alterations in the gut microbiota have been linked with compromised gut barrier function, immune dysregulation, and systemic spread of the gut microbiota, as discussed by Jones et al. [52] and highlighted above.However, SCFA production as not one of the metabolic pathways associated with growth or growth velocity in the only metagenomic study in children living in rural Zambia [51].
Fig 8 shows the Spearman correlation coefficients for the taxa that significantly correlated to the fecal SCFA concentrations.Most taxa correlated to propionate, either positively or negatively, even if those taxa are not considered to be directly involved in propionate production, such as Faecalibacterium, which produces butyrate.None of the other studies on stunting or growth have looked at SCFA in feces.
The strength of this study is the inclusion of 100 children in each group.Most previous studies have only looked at a limited numbers of stunted children [47][48][49][50][51].However, the study also has limitations.All taxa currently at most correlate to being stunted; no cause and effect relationship has been established.Also, as is clear from our own two studies in Indonesia (on Java and ENT), as well as the study by Rinanda et al. [47] on Sumatra, the composition of the gut microbiota, even within one country, is geography specific, which is corroborated by the studies in Zimbabwe, India and Peru [47][48][49][50][51]. Hence no general conclusions can be drawn.On top of that, differences in library preparation, sequencing primers, PCR cycles or bioinformatics tools may lead to different results [53].Although for our study on Java and ENT different libraries were prepared (in fact even within those studies more than 1 library was necessary due to the number of samples), all other tools were identical.Based on our positive control and mock control there were no differences between the sequence runs for those two geographical location, and no other technical difference were present (primers, number of PCRcycles, etc.).As such, the difference in composition between the gut microbiota on Java and ENT is believed to be primarily due to geographical location, and influenced by differences in diet, hygienic status and environmental conditions.

Conclusion
In conclusion, the relationship between gut microbiota and stunting is complex, and is influenced by many confounding factors.In the current study three taxa were consistently different  between stunted children and non-stunted children at both sampling sites in East Nusa Tenggara.However, these were different taxa compared to our study on Java, even though this was in the same country, partly as a consequence of the different composition of the samples from ENT altogether.Differences in diet, cultural habits, but certainly also infrastructure facilities such as clean water and sanitation are expected to contribute to the differences between our studies on Java and ENT, as we show here that these sanitary parameters are different between the two groups, and correlate to several taxa.Based on recent results from Robertson at el. [51] on whole metagenome sequencing, it may not be sufficient to indicate differences in taxa, but important biochemical pathways may have to be elucidated, including those for SCFA production, which link to barrier function and immune modulation.Further, it remains to be seen whether the taxa that positively correlated with zlen are involved in stunting, and whether increasing the ones that are negatively correlated with zlen might correct some of the growth retardation.In the meantime, we are performing a food intervention on these children to investigate this.

Fig 1 .
Fig 1. CONSORT flow diagram of the cross-sectional study in East Nusa Tenggara, Indonesia.https://doi.org/10.1371/journal.pone.0299349.g001 Fig 3A shows the principal coordinate analysis (PCoA) of the Bray-Curtis dissimilarity between the samples, colored by being stunted or not.The PCoA plot for the other indices are shown in S1 Fig.

Fig 3 .
Fig 3. Principal coordinate analysis plot of Bray Curtis dissimilarity for A) non-stunted (red) and stunted (blue) children of the current study in East Nusa Tenggara, B) the non-stunted (red) and stunted (blue) children of the current study in East Nusa Tenggara (full circles) and our previous study on Java (open circles).https://doi.org/10.1371/journal.pone.0299349.g003

Fig
Fig 6A shows the taxa that are significantly (q < 0.05) correlated (Spearman) with height, weight and bmi, and Fig 6B with the corresponding z-scores.The cells are colored by the correlation coefficient (rho).As expected from the significantly higher RA of Faecalibacterium, Lachnoclostridium and Veillonella, these are negatively correlated with height and zlen.As expected there is significant overlap between the taxa in Fig 6A and the metrics derived from the anthropometric indices in Fig 6B, although there are a few differences.Most taxa are correlated to bmi and zbmi (and the similar zwfl).Only an uncharacterized genus in the Ruminococcacea family correlates negatively to weight and zwei.Besides Faecalibacterium, Lachnoclostridium and Veillonella, there are 5 (Fig 6A) and 9 (Fig 6B) other taxa that are negatively correlated to height and zlen, respectively, while there are 14 (Fig 6A) and 13 (Fig 6B) taxa that are positively correlated to height and zlen, respectively.

Fig 8 .
Fig 8. Heatmap of the taxa at genus level that were different when correlated by Spearman correlation analysis (q < 0.05) to the fecal concentrations of the SCFAs acetate, propionate and butyrate.Cells are colored by rho (correlation coefficient): Blue: Positive correlation; red: Negative correlation.https://doi.org/10.1371/journal.pone.0299349.g008 Food Innovation (HEFI) of Maastricht University -campus Venlo, which has been made possible with the support of the Dutch Province of Limburg; d) Ilia Popov was funded by the Russian Federation through the scholarship of the President of the Russian (Order of the Ministry of Science and Higher Education of the Russian Federation No. 825 of 29 August 2022).The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Table 2 . Sanitary parameters that were significantly different between the groups †. stunted normal p-value
'around', or 4: Private, inside), handwashing after toilet visit (1: No or 2: Yes), handwashing with flowing water (1: No or 2: Yes), handwashing with soap and/or (flowing) water (1: Neither, 2: Flowing water, 3: Water and soap, or 4: Flowing water and soap).Arbitrary values of the non-stunted group (normal) were set to 1, and values for the stunted group were subsequently normalized to the normal group.https://doi.org/10.1371/journal.pone.0299349.t002